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Abstract 

A dynamical system model for tumor - immune system interaction together 
with a method to mimic radiation therapy are proposed. A large population of 
virtual patients is simulated following an ideal radiation treatment. A characteristic 
parameter, the Immune System - Tumor Efficiency Rate (ISTER), is introduced. 
ISTER dependence of treatment success and other features is studied. Statistical 
results allow us to give a patient classification scheme. Radiotherapy treatment 
biological effective dose (BED) is thus optimized based on the patient physical 
condition, following the ALARA (As Low As Reasonably Achievable) criterion. 

1 Introduction 

Some approaches to cancer growth and behavior have been made in the past years. 
Recent techniques try to use a population dynamics model JT] [2] [3] |4) to mathemati- 
cally describe the tumour behavior and its interaction with the immune system. Some 
of these works explain tumour behavior under clinical treatments like cytokines Q or 
radiovirotherapy |6j and properly explain the qualitative behaviors of several tumours. 
Even though great efforts had been made to describe cancer radiotherapy treatments 
13, they are but vaguely linked to clinical observations and their large number of vari- 
ables and coefficients make their results hardly transposable to a clinical context. 

Radiotherapy and surgery are the most effective treatments for cancer, and even 
while surgery has a longer tradition, radiotherapy is replacing surgery for the control 
of many tumours 1 7 j . Those treatments follow strict protocols that often apply a fixed 
physical radiation dose, hardly taking into account the kind of tumour or the patient 
immunological condition. Thus, a radiotherapy protocol might result in a very low 
success probability for some patients starting their treatments with a weakened im- 
mune system. In practice such a treatment will be interrupted if the patient physical 
condition worsens, although the patient will have already received inappropriate doses 
of radiation. 

Due to its importance and looking for an applicable method, we intend to model 
a radiotherapy treatment making the simplest possible assumptions. Furthermore we 
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will introduced the Immune System Tumour Efficiency Rate parameter (ISTER) as a 
measure of the patient immune system strength to fight back cancer. This parameter 
allow us to make a patient classification and find the success probability of each patient 
group following a radiotherapy treatment protocol. Finally, we will use these results 
to assess the optimized biological effective dose (BED) or tissue effect (E) based on a 
given patient physical condition. 

2 Model 

In order to describe the tumour evolution, we propose a Lotka-Volterra like model 
based on some assumptions. Tumour cells growth X (as usual, a dot over a quantity 
represents its time derivative) depends on the current tumour population as aX and its 
mass-law interaction with lymphocytes, —bXY. Lymphocytes population grows due to 
tumour-immune system interaction, dXY, and falls in time exponentially, —fY, due to 
natural cell death. Tumour secretes interleukin which produces an immune depression 
effect |8]|9], an d we will make the simplest assumption supposing it proportional to 
the tumour cell number,— kX. The tumour is localized and there is a constant flow, u, 
of lymphocytes from the immune system into this region. 

So, we will model tumour-immune system interaction using the known equations 

0: 

X = aX-bXY 

Y = dXY -fY -kX + u { > 

The effects of radiation over any tissue are generally classified in three phases 
[7). Physical phase, when radiation ionizes atoms. Chemical phase, when ionized 
molecules interact with other biological components of the cell. And finally, biolog- 
ical phase, where the damage is fixed, and unrepairable cells are signaled to die by 
apoptosis. 

Carcinogenesis and other malignant effects, that escape cellular control, can appear 
as late effects of the biological phase and, to avoid them, radiation doses need to be op- 
timized. This means that higher doses that could reduce the long term overall survival 
of patients ifTTI . must be avoided whenever possible. 

We will collect all these heterogeneous effects, according to their time scale, in 
two groups: short and long term effects. Short term effects occur at very small time 
scales compared with the time scales on which our model runs (those times for which 
changes in the X and Y variables become appreciable), and so only long term effects 
will be taken into account in our evolution equations. Then, we are going to assume 
that lymphocytes die or loose their ability to attack tumour cells immediately, and that 
radiation dose is concentrated at an infinitesimal instant of time. At that very moment, 
long term effects start to take place, whereas short term effects instantaneously modify 
the state of the system. 

Thus, we also assume that when a radiation dose is applied at a given instant T„, 
it induces a fraction B t of the tumour cells to lose their reproductive endowment and 
to die exponentially. The fraction S t of tumour cells not affected by radiation can be 
computed by the linear-quadratic (LQ) model Bl ITOl . 
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S, = 1 -B t = exp[-£] = exp[-aA - pA 2 ] (2) 

where E is known as the tissue effect, a and p are Type A and B damage coefficients 
171 , and A is the physical radiation dose expressed in Gy, as usual in clinical contexts. 
Furthermore, a fraction B/ of lymphocytes is also killed by radiation, in a manner sim- 
ilar to © although having different a and p coefficients. 

To include long term processes in Eqs. (Q]), we write a new equation for non- 
proliferating tumour cells Z (5), taking into account that lymphocyte population is also 
stimulated, as pZY, due to interaction with these cells. The number of non-proliferating 
tumour cells decays exponentially as — rZ due to the death of damaged cells, and also 
as —qZY due to the interaction with lymphocytes. Then we arrive to the system 

X = aX-bXY-B,(T)X 

Y = dXY + pZY -fY- k(X +Z)+u- B,(T)Y (3) 
Z = B,(T)X-rZ-qZY 

where B t (T) = B,Y.8(T - T„) and B,(T) = B/£8(r - T n ). T n are the time instants 
when radiation doses are applied and 8(T — T„) denotes Dirac's delta centered at T n . We 
have supposed that lymphocytes interact in different ways with X and Z cells, although 
both kind of tumour cells cause the same depression over the immune system. 

Equations ([3]) can be expressed in a dimensionless form taking the tumour dupli- 
cation time x c = 1/a (in absence of external influences) as the characteristic time, so 
we introduce the dimensionless time x = T /x e . Through the substitutions X = ax/d, 
Y = ay/b, Z = az/d, we obtain the dimensionless system: 

i = x — xy — J t (x)x 

y = xy + ezy - Xy - k(x + z) + o - y/(x)y (4) 
z = y t (i)x - pz - y\zy 

with y/ = B{/a, y t = B t /a, e = p/d, A, = / /a, K = kb/ ad, a = ub/a 2 , p = ra/d and 
T| = qa 2 /db. 

A linear stability analysis of the system dU shows that tumour will vanish to Lq = 
(0;o/X;0) if a/A, > 1 and will remain controlled around L\ = ((A, — ct)/(1 — k);1;0) 
ifK<a/A< 1 orK> 1 13- If the system is Lo-stable and initial tumour size is small 
enough, then the radiation treatment is unnecessary, whereas if tumour size is large 
enough, then the treatment will take it closer to Lq. 

The L\ controlled growth state will be reached only if both parameters fulfill the 
same condition, in other words, if a/A and K are both greater or smaller than unity 
at the same time. Any other condition makes L\ < 0, and even when the stable point 
mathematically exists, it can not be approximated from realistic initial conditions (that 
should remain positive along the simulation time). For those patients with K > 1 and 
a/A, < 1, the main effects of the tumour will be the depression of immune system, 
they will present a low Karnofsky performance scale ITD and will not fulfill physical 
conditions to be subject under treatment. 

However for a/A < K < 1, tumour will grow exponentially and radiotherapy goal 
will be to bring it close enough to Lq so that immune system can get rid of the tumour. 
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Figure Q] shows stable and unstable regions of Eqs. and highlights region III on 
which this work will focus. 
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Figure 1: Phase diagram of equations @] The focus of this work falls inside shadowed 
region. 

The chosen characteristic time and the dimensionless parameters allow us to give a 
very intuitive interpretation of the critical parameters of Eqs. We can see a/A, as 
the efficiency of immune system over tumour growth and K as the "deficiency" of the 
immune system due to tumour growth. 

It is also easy to see that radiation treatments do not change the stability conditions 
of our system, given that radiotherapy does not change tumour or lymphocytes growth 
rate, but can drive the number of both kind of cells to very small values. Although 
Eqs. (01) allow for infinitesimal x values, in real systems when the number of tumour 
cells becomes small enough, immune system may kill them J7). However, in other 
cases when a few tumour cells survive, they can cause tumour regrowth. It is known 
that this behavior is almost independent on tumour size and as an estimation we will 
assume that the closer Lq is (in terms of the phase space of figure Q]) to the line where 
it becomes an stable point, the higher will be the probability of tumour elimination by 
the immune system. Thus, when x becomes small enough we will take 



as the probability of tumour regression. 

Similarly, whenever lymphocyte population becomes zero we will assume a general 
failure of immune system. This situation may also occur in some cases where the 




(5) 
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tumour is removed but the immune system reaches such an extreme low concentration 
of lymphocytes that, consequently, the patient dies. 

3 Simulation 

We can mimic different radiation treatments with Eqs. (HJi to simulate tumour evolu- 
tion. To follow radiotherapy treatment in a realistic way, we apply a radiation session 
every workday and none in weekends. All treatments lfl2l[T3l begin the tenth day, take 
6 weeks of radiotherapy and patients are under observation until 6 months after the end 
of radiotherapy sessions. We generate several virtual patients under treatment taking 
different values for the parameter values in Eqs. |4]and use a four step Runge-Kutta 
method |[T4l to integrate them. 

To reproduce tumour evolution resembling that of a clinical case, we need to cal- 
culate the correct values of the coefficients appearing in Eqs. (Q]). Estimation of these 
coefficients was made in 1131 , and following a similar procedure it would not be too 
hard for clinical professionals to estimate their values. Figure [2] shows treatment evo- 
lution for typical values of the coefficients and under different doses of radiation. We 
can see how the number of tumour cells capable of mitosis quickly decreases with the 
radiation therapy. For long enough times, if regression behavior is not accomplished, 
the tumour regrows quickly. 

In order to accomplish an statistical study of the dependence of treatment success 
on the dosage, and due to the wide range of possible parameter values in Eqs. |@}, their 
values are drawn randomly from a log-normal distribution, to avoid negative values, 
but keeping the efficiency of immune system {o/X) always smaller than 1. Survival 
factors 131|2l are also taken as random values within the interval shown in table Q] As 
initial conditions we have supposed, for simplicity, that the number of tumour cells 
is higher than the number of lymphocytes and that both populations are distributed as 
normal random numbers, with parameters shown in table 12 We have also tested other 
distributions for the initial conditions as well as for coefficient values, to verify that the 
choice does not affect the qualitative nature of our results. 



Parameter 


Minimum 


Maximum 


X 


10° 


10 J 


a 


io- 1 


IO 5 


K 


io- 2 


10 4 


St 


0.5 


0.9 


s, 


0.1 


0.4 



Table 1: Dimensionless parameter values of equations [4] taken from IPT51 l4l [Toll . 

At this point we can proceed to make statistical predictions by generating a popula- 
tion of "virtual patients" (characterised by their immune system and tumour parameter 
values) and simulating their treatment evolutions. Tables [3] and [2] show parameter val- 
ues used to generate virtual patients. 
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S,= 0.8 



10 20 30 40 50 

T (days) 

Figure 2: Tumour evolution under radiotherapy treatment for two different tumour 
survival factors. A, = 5.0, O = 4.0, K = 0.7, S, = 0.2 



Coefficient 


Mean 


Standard deviation 


s, 


0.6 


0.1 


Si 


0.18 


0.06 


x 


1.0 


0.1 


yo 


0.5 


0.1 



Table 2: Statistical survival factors and initial conditions for tumour cells and T- 
lymphocytes 



Parameter 


log. Mean 


log. Standard Deviation 


X 


5.0 


0.5 


a 


2.5 


0.5 


K 


0.8 


0.2 



Table 3: Statistical parameter values. 
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4 Results and clinical interpretation 



We have created a database consisting of over 3 x 10 5 virtual patients. We have calcu- 
lated the probability of treatment success (P s ) as the fraction of patients without tumour 
at the end of treatment. We have represented this probability P s as a function of tissue 
effect, E, (see equation (fJI) and efficiency of immune system (a/A,) or ISTER. In Fig 
[3] a color map of P s versus E and a/A is represented. This allows us to classify patients 
based on their c/A-value and to assess those patients to whom, having an extremely 
low success probability, the application of high radiation doses would render useless. 
Radiotherapy is not the appropriate treatment for those patients, although it could be 
used as a palliative, if a good balance between drawbacks and advantages is presumed 
for a specific patient. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



a/A, 

Figure 3: P s as a function of E and a/A. E- curve corresponds to the minimum values 
of tissue effect for which P s > 0. E + curve represents the values of tissue effect for 
which P s reaches its maximum value, given a fixed a/A. Marks E\ = 0.84 and Ei = 
0.24 represents the tissue effect in the explained examples. 

We can see that, for a given value of a/A, two significant values of E can be defined: 
below which P s is very small (less than 1%), and E + , above which P s is almost 
constant (with less than 1% of change). Results can be fitted to the expression, 

P s = P s (a/X,E) (6) 

and the significant values of E computed as functions of a/A. These two threshold 
values (£_ and E + ) divide the phase space (a/A,£) into three regions as shown in Fig 
[3] The success probability is negligible in region /, below while it almost attains 
its maximum value above E + , in region III . However, on the intermediate region 77, as 
E grows, P s increases faster towards its maximum value (above the E + curve). 
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The coefficients a and p are generally hard to find, but not impossible, and several 
values of the ratio oc/p are reported in the literature Q. However this is not enough 
for the clinical application and at least one of them must be found (as explained also in 
|7|) to proceed. Luckily, to characterize patients, we just need, among all coefficients 
involved in Eqs. (0), to know the ratio u/f, the effective amount of lymphocytes in the 
absence of tumour effects or inmunodepression, and b/a, a measure of the effective- 
ness of lymphocytes over tumour growth. Clinical professionals must determine the 
inquiries and tests needed to find a patient's ISTER. 

To illustrate a possible clinical application of this result, we are going to suppose 
two virtual patients with the same ISTER = 0.5 and different tumour sensitivities. We 
will assume a sensitive tumour 11611 with a = 0.3 Gy 1 and oc/p = 5 Gy, and a more 
resistant tumour with a = 0.1 Gy 1 and oc/p = 10 Gy, so the tumour resistance to 
radiation is quite different for each case. In both cases the tissue effects, with an usual 
treatment, are represented in Fig [3] like E\ and E% respectively. 

First, let us consider the case of a sensitive tumour lfl6l with a = 0.3 Gy -1 and 
Oc/P = 5 Gy. If we apply the typical fractionated radiotherapy used in our calculations, 
then the biological effective dose {BED = E /a), for each radiotherapy session of 2 Gy, 
will be 2.8 Gy. However, this high dose value does not really increase the success 
probability. A patient with an immune system efficiency of ISTER = 0.5, has his 
maximum healing probability for a BED value around 1.5 Gy in each radiotherapy 
session. Then, we must apply a physical radiation of 1 .2 Gy in each radiation session 
and thus, avoid an useless amount of 24 Gy to be applied in the whole treatment. 

However, in the case of a tumour having a higher resistance to radiation, e.g. with 
a = 0. 1 Gy~' and oc/p = 10 Gy, the same patient with an ISTER = 0.5, attains a max- 
imum success probability for a BED = 4.5 Gy, and needs a physical dose of 3.4 Gy to 
be applied in each session. Besides, the minimum BED value is 2.2 Gy, correspond- 
ing to 1.85 Gy of physical radiation per session. Table|4]shows the optimal calculated 
values of physical radiation dose for the two examples of tumour with a ISTER = 0.5. 

The oncologist, should decide the amount of radiation to apply, by evaluating Eq. 
(O to know the treatment success probability, and taking into account any other clinical 
factors implied. 



Kind of tumour 


BED for 

A = 2 Gy 


Optimal E 


Optimal 
BED 


Optimal A 
per session 


Sensitive tumour 
oc = 0.3Gy 
oc/p=5Gy 


2.8 Gy 


0.45 


1.5 Gy 


1.2 Gy 


Less sensitive tumoui 
a = 0.1 Gy 
a/p = 10 Gy 


2.4 Gy 


0.45 


4.5 Gy 


3.4 Gy 



Table 4: Optimal values for two different tumours with a rate ISTER = 0.5. 



The presented results match with those reported in ifTTI . that show that the long term 
survivance of patients is not better at higher doses of radiation. On the contrary, the 
higher number of long term survival patients is reached at intermediate doses (between 
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2 or 3 Gy), even with a smaller total amount of radiation. 

5 Conclusions 

The proposed method, allows us to find the success probability of a fractionated radio- 
therapy treatment, using the patient ISTER parameter, as a new oncological index, and 
the survival fraction S t of tumour cells, even if other parameters involved are unknown. 
This calculation provides a way to classify patients, based on their ISTER value, and 
to approach to the optimum treatment. 

The radiotherapy treatment must be designed for each patient taking into account 
his/her immunological characteristics (ISTER) relative to the tumour. Tissue effect has 
to be tuned to be larger than otherwise no success will be achieved, but needs not to 
be larger than E+, because no improvement will be obtained for larger radiation doses. 
Thus, in accordance with the ALARA (As Low As Reasonably Achievable) principle 
ifTTl . the physical radiation doses should be adjusted to bring E as close as possible 
to E+ but without out-ranging it. This optimization process could be performed once 
the clinical professionals find a way to evaluate the ISTER index experimentally for a 
given patient. On other hand, the values of of a and p (in Eq. |2} are known or feasible 
to find for many kinds of tumour. 
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A Fitting result data to an analytical function 



We used a Levenberg-Marquardt lfl4l method to fit the result data, showed in figure[3] 
to the analytical function, 



fi(qA,£) 
a/X 




(7) 



for each value of computed a/A,. This expression gives us a family of functions related 
to each other through coefficients 0, 0, \\t and cp. These coefficients are functions of 
only a/a and can be easily fitted using the same numerical method. 

We have found the following numerical expressions for these coefficients, 

9 = 0.950271 * (l-exp(-4.66627£- 0.24319)) 
(j) = -0.935012 + exp(-4.71719f - 0.289458) 

<p = 0.0450091 f + 0.091267 () 
\|/ = 0.0159581^ -0.141425 

Merging all this expressions, it is possible analyze the behaviour of the success 
probability P s (o/X,E). 
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